This file and all other referenced in the code can be found at the repository: https://github.com/Andros-Spica/Quaternary_Angourakis-et-al-2021
Load packages:
require(reshape2)
## Loading required package: reshape2
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
require(patchwork)
## Loading required package: patchwork
## Warning: package 'patchwork' was built under R version 3.6.3
require(colorspace)
## Loading required package: colorspace
## Warning: package 'colorspace' was built under R version 3.6.3
require(grid)
## Loading required package: grid
require(gridExtra)
## Loading required package: gridExtra
## Warning: package 'gridExtra' was built under R version 3.6.3
require(png)
## Loading required package: png
Define lookup functions to access data between tables using identifiers or ‘keys’:
getRelatedVariableWithKey <- function(keyInCallingTable, keyInCalledTable, variableInCalledTable)
{
return(
sapply(keyInCallingTable, function(x) variableInCalledTable[match(x, keyInCalledTable)])
)
}
getRelatedVariableWithTwoKeys <- function(firstKeyInCallingTable, secondKeyInCallingTable,
firstKeyInCalledTable, secondKeyInCalledTable, variableInCalledTable)
{
temp <- 1:length(firstKeyInCallingTable)
for (i in 1:length(temp))
{
temp[i] <- variableInCalledTable[firstKeyInCalledTable == firstKeyInCallingTable[i] &
secondKeyInCalledTable == secondKeyInCallingTable[i]]
}
return(temp)
}
Define the format of image files to be generated (only “png” and “eps”):
listOfImageFormats <- c("eps", "png")
Load crop table:
cropTable <- read.csv("models/cropsTable.csv", skip = 4)
Clean empty columns (generated with automatic names starting with X):
cropTable <- cropTable[,names(cropTable)[grep("^(?!X).*$", names(cropTable), perl = TRUE)]]
Reorder crops for a more logical display:
WARNING: doing this manually makes the script sensitive to the value of crop-selection in NetLogo. Thus, this chunk must be updated if crop-selection is changed before running experiments.
cropTable$crop <- factor(cropTable$crop, levels = c("proso millet", "pearl millet", "rice", "barley", "wheat 1", "wheat 2"))
Get vector of colours to represent crops:
cropColours <- qualitative_hcl(nlevels(cropTable$crop))
# alternative using only base R graphics
#cropColours <- rainbow(nlevels(cropTable$crop), s = 0.8, v = 0.8, end = 0.85)
knitr::kable(cropTable)
| crop | species | cultivar | T_sum | HI | I_50A | I_50B | T_base | T_opt | RUE | I_50maxH | I_50maxW | T_heat | T_ext | S_CO2 | S_water | sugSowingDay | sugHarvestingDay | root.zone.depth..mm. |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| wheat 1 | Triticum aestivum | Yecora Rojo | 2200 | 0.36 | 480.0 | 200.00 | 0 | 15 | 1.24 | 100 | 25 | 34 | 45 | 0.08 | 0.40 | 330 | 105 | 1000 |
| wheat 2 | Triticum aestivum | Batten | 2150 | 0.34 | 280.0 | 50.00 | 0 | 15 | 1.24 | 100 | 25 | 34 | 45 | 0.08 | 0.40 | 330 | 105 | 1000 |
| rice | Oryza sativa | IR72 | 2300 | 0.47 | 850.0 | 200.00 | 9 | 26 | 1.24 | 100 | 10 | 34 | 50 | 0.08 | 1.00 | 170 | 320 | 400 |
| barley | Hordeum vulgare | hypothetical; based on DSSAT v4.7 | 1762 | 0.42 | 350.0 | 170.00 | 0 | 15 | 1.24 | 100 | 20 | 34 | 45 | 0.08 | 0.40 | 330 | 100 | 1000 |
| pearl millet | Pennisetum glaucum | hypothetical; based on DSSAT v4.7 - MLCER047.CUL - Middle Variety | 1220 | 0.25 | 245.0 | 120.00 | 10 | 33 | 1.90 | 100 | 5 | 35 | 47 | 0.08 | 0.05 | 200 | 320 | 1000 |
| proso millet | Panicum miliaceum L. | hypothetical | 1328 | 0.29 | 157.3 | 96.75 | 0 | 18 | 3.00 | 100 | 130 | 34 | 45 | 0.08 | 0.05 | 107 | 220 | 1000 |
Experiment 1 simulates the growth and yield of crops included in ‘cropTable’ using the NetLogo implementation of SIMPLE crop model (Zhao et al. 2019) adapted as a module of the Indus Village model under the default parameter configuration, which is aimed to approximate the conditions in Haryana, NW India.
Load data:
yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")
Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData$yield[is.na(yieldData$meanARID_grow)] <- NA
Force the same order of crops used by cropTable:
yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))
Calculate the range of ARID and productivity (yield):
minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)
Load daily simulated weather and ARID data:
weatherData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_weather_type-of-experiment=user-defined_experiment-name=exp1_initRandomSeed=0.csv")
As the reference real-world data, we use the data downloaded at NASA´s POWER access viewer (power.larc.nasa.gov/data-access-viewer/) selecting the user community ‘Agroclimatology’ between 01/01/1984 and 31/12/2007. The exact location used as a reference is:
We selected the ICASA Format’s parameters:
Load Soil Water (or Soil Water Balance) model code to calculate ARID given real-world data:
source("FigWeather_input/soilWaterModel/estimateETr.R")
source("FigWeather_input/soilWaterModel/watbal.model.R")
Load and format the real-world dataset so that they can be used in the Soil Water model:
weatherData_rakhigarhi <- watbal.weather.file("FigWeather_input/POWER_SinglePoint_Daily_19840101_20071231_029d17N_076d70E_5b401917.csv",
year = 1984:2007)
Estimate reference evapotranspiration using FAO recomendations (Penman-Monteith method):
weatherData_rakhigarhi$ETr <- estimateETr(weatherData_rakhigarhi$I,
weatherData_rakhigarhi$T2M,
weatherData_rakhigarhi$Tmax,
weatherData_rakhigarhi$Tmin,
weatherData_rakhigarhi$T2MDEW,
weatherData_rakhigarhi$WS2M)
Calculate ARID in real-world data:
Use Soil Water model to calculate soil water content proportion or WATp (mm mm-1 day-1) and ARID coefficient.
The first step is to set all soil parameters, here assumed to be constant (the same values used in experiment 1 simulations):
watbalParameters <- watbal.define.param()
Run the model inputing parameters and dataset:
watbalOutput_rakhigarhi <- watbal.model(watbalParameters,
weatherData_rakhigarhi,
typeOfParValue = "nominal")
watbal.modeliterates for every day in the dataset calculating soil water content, absolute and proportion (WAT, WATp), and the respective ARID coefficient.
Show parameter values used here in Soil Water model (both in watbal.model and in the SIMPLE crop model). In experiemnt 1, they should be the same:
knitr::kable(
cbind(
parameters = c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF"),
SIMPLE = as.vector(unlist(yieldData[1, c("elevation", "albedo", "CN", "DC", "z", "FC", "WHC", "WP", "MUF")])),
watbal = c(200, # these values of elevation and albedo are the default in estimateETr
0.23,
as.vector(watbalParameters[1, c("CN", "DC", "z", "FC", "WHC", "WP", "MUF")]))
)
)
| parameters | SIMPLE | watbal |
|---|---|---|
| elevation | 200 | 200 |
| albedo | 0.23 | 0.23 |
| CN | 65 | 65 |
| DC | 0.55 | 0.55 |
| z | 400 | 400 |
| FC | 0.21 | 0.21 |
| WHC | 0.15 | 0.15 |
| WP | 0.06 | 0.06 |
| MUF | 0.096 | 0.096 |
Create data frames containing the daily summary statistics of simulations and real-world data:
weatherSummaryNames <- c("dayOfYear",
"solarRadiation.mean",
"solarRadiation.sd",
"solarRadiation.max",
"solarRadiation.min",
"solarRadiation.error",
"temperature.mean",
"temperature.sd",
"temperature.max",
"temperature.min",
"temperature.error",
"maxTemperature.mean",
"maxTemperature.max",
"maxTemperature.min",
"maxTemperature.error",
"minTemperature.mean",
"minTemperature.max",
"minTemperature.min",
"minTemperature.error",
"temperature.lowerDeviation",
"temperature.lowerDeviation.error",
"temperature.upperDeviation",
"temperature.upperDeviation.error",
"precipitation.mean",
"precipitation.max",
"precipitation.min",
"precipitation.error",
"ARID.mean",
"ARID.max",
"ARID.min",
"ARID.error")
weatherSummary <- vector("list", length(weatherSummaryNames))
names(weatherSummary) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential
for (day in 1:365)
{
weatherSummary$dayOfYear <- c(weatherSummary$dayOfYear, day)
tempData_thisDay <- weatherData[weatherData$currentDayOfYear == day,]
# solar radiation
weatherSummary$solarRadiation.mean <- c(
weatherSummary$solarRadiation.mean,
mean(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.sd <- c(
weatherSummary$solarRadiation.sd,
sd(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.max <- c(
weatherSummary$solarRadiation.max,
max(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.min <- c(
weatherSummary$solarRadiation.min,
min(tempData_thisDay$solarRadiation, na.rm = TRUE))
weatherSummary$solarRadiation.error <- c(
weatherSummary$solarRadiation.error,
qt(0.975,
length(tempData_thisDay$solarRadiation) - 1) *
sd(tempData_thisDay$solarRadiation, na.rm = TRUE) /
sqrt(length(tempData_thisDay$solarRadiation)))
# temperature
## daily mean
weatherSummary$temperature.mean <- c(
weatherSummary$temperature.mean,
mean(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.sd <- c(
weatherSummary$temperature.sd,
sd(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.max <- c(
weatherSummary$temperature.max,
max(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.min <- c(
weatherSummary$temperature.min,
min(tempData_thisDay$temperature, na.rm = TRUE))
weatherSummary$temperature.error <- c(
weatherSummary$temperature.error,
qt(0.975,
length(tempData_thisDay$temperature) - 1) *
sd(tempData_thisDay$temperature, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature)))
## daily max
weatherSummary$maxTemperature.mean <- c(
weatherSummary$maxTemperature.mean,
mean(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.max <- c(
weatherSummary$maxTemperature.max,
max(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.min <- c(
weatherSummary$maxTemperature.min,
min(tempData_thisDay$temperature_max, na.rm = TRUE))
weatherSummary$maxTemperature.error <- c(
weatherSummary$maxTemperature.error,
qt(0.975,
length(tempData_thisDay$temperature_max) - 1) *
sd(tempData_thisDay$temperature_max, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_max)))
## daily min
weatherSummary$minTemperature.mean <- c(
weatherSummary$minTemperature.mean,
mean(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.max <- c(
weatherSummary$minTemperature.max,
max(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.min <- c(
weatherSummary$minTemperature.min,
min(tempData_thisDay$temperature_min, na.rm = TRUE))
weatherSummary$minTemperature.error <- c(
weatherSummary$minTemperature.error,
qt(0.975,
length(tempData_thisDay$temperature_min) - 1) *
sd(tempData_thisDay$temperature_min, na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_min)))
## daily upper and lower deviation
weatherSummary$temperature.lowerDeviation <- c(
weatherSummary$temperature.lowerDeviation,
mean(tempData_thisDay$temperature -
tempData_thisDay$temperature_min))
weatherSummary$temperature.lowerDeviation.error <- c(
weatherSummary$temperature.lowerDeviation.error,
qt(0.975,
length(tempData_thisDay$temperature_min) - 1) *
sd(tempData_thisDay$temperature -
tempData_thisDay$temperature_min,
na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_min)))
weatherSummary$temperature.upperDeviation <- c(
weatherSummary$temperature.upperDeviation,
mean(tempData_thisDay$temperature_max -
tempData_thisDay$temperature))
weatherSummary$temperature.upperDeviation.error <- c(
weatherSummary$temperature.upperDeviation.error,
qt(0.975,
length(tempData_thisDay$temperature_max) - 1) *
sd(tempData_thisDay$temperature_max -
tempData_thisDay$temperature,
na.rm = TRUE) /
sqrt(length(tempData_thisDay$temperature_max)))
# precipitation
weatherSummary$precipitation.mean <- c(
weatherSummary$precipitation.mean,
mean(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.max <- c(
weatherSummary$precipitation.max,
max(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.min <- c(
weatherSummary$precipitation.min,
min(tempData_thisDay$RAIN, na.rm = TRUE))
weatherSummary$precipitation.error <- c(
weatherSummary$precipitation.error,
qt(0.975,
length(tempData_thisDay$RAIN) - 1) *
sd(tempData_thisDay$RAIN, na.rm = TRUE) /
sqrt(length(tempData_thisDay$RAIN)))
# ARID
weatherSummary$ARID.mean <- c(
weatherSummary$ARID.mean,
mean(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.max <- c(
weatherSummary$ARID.max,
max(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.min <- c(
weatherSummary$ARID.min,
min(tempData_thisDay$ARID, na.rm = TRUE))
weatherSummary$ARID.error <- c(
weatherSummary$ARID.error,
qt(0.975,
length(tempData_thisDay$ARID) - 1) *
sd(tempData_thisDay$ARID, na.rm = TRUE) /
sqrt(length(tempData_thisDay$ARID)))
}
rm(day, tempData_thisDay)
weatherSummary <- data.frame(weatherSummary)
weatherSummary_rakhigarhi <- vector("list", length(weatherSummaryNames))
names(weatherSummary_rakhigarhi) <- weatherSummaryNames
# OBS: the lines above produce an ERROR related to names that is inconsequential
for (day in 1:366)
{
weatherSummary_rakhigarhi$dayOfYear <- c(weatherSummary_rakhigarhi$dayOfYear, day)
tempData <- weatherData_rakhigarhi[weatherData_rakhigarhi$day == day,]
tempWatbalData <- watbalOutput_rakhigarhi[watbalOutput_rakhigarhi$day == day,]
# solar radiation
weatherSummary_rakhigarhi$solarRadiation.mean <- c(
weatherSummary_rakhigarhi$solarRadiation.mean,
mean(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.sd <- c(
weatherSummary_rakhigarhi$solarRadiation.sd,
sd(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.max <- c(
weatherSummary_rakhigarhi$solarRadiation.max,
max(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.min <- c(
weatherSummary_rakhigarhi$solarRadiation.min,
min(tempData$I, na.rm = TRUE))
weatherSummary_rakhigarhi$solarRadiation.error <- c(
weatherSummary_rakhigarhi$solarRadiation.error,
qt(0.975, length(tempData$I) - 1) *
sd(tempData$I, na.rm = TRUE) /
sqrt(length(tempData$I)))
# temperature
## daily mean
weatherSummary_rakhigarhi$temperature.mean <- c(
weatherSummary_rakhigarhi$temperature.mean,
mean(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.sd <- c(
weatherSummary_rakhigarhi$temperature.sd,
sd(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.max <- c(
weatherSummary_rakhigarhi$temperature.max,
max(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.min <- c(
weatherSummary_rakhigarhi$temperature.min,
min(tempData$T2M, na.rm = TRUE))
weatherSummary_rakhigarhi$temperature.error <- c(
weatherSummary_rakhigarhi$temperature.error,
qt(0.975, length(tempData$T2M) - 1) *
sd(tempData$T2M, na.rm = TRUE) /
sqrt(length(tempData$T2M)))
## daily max
weatherSummary_rakhigarhi$maxTemperature.mean <- c(
weatherSummary_rakhigarhi$maxTemperature.mean,
mean(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.max <- c(
weatherSummary_rakhigarhi$maxTemperature.max,
max(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.min <- c(
weatherSummary_rakhigarhi$maxTemperature.min,
min(tempData$Tmax, na.rm = TRUE))
weatherSummary_rakhigarhi$maxTemperature.error <- c(
weatherSummary_rakhigarhi$maxTemperature.error,
qt(0.975, length(tempData$Tmax) - 1) *
sd(tempData$Tmax, na.rm = TRUE) /
sqrt(length(tempData$Tmax)))
## daily min
weatherSummary_rakhigarhi$minTemperature.mean <- c(
weatherSummary_rakhigarhi$minTemperature.mean,
mean(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.max <- c(
weatherSummary_rakhigarhi$minTemperature.max,
max(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.min <- c(
weatherSummary_rakhigarhi$minTemperature.min,
min(tempData$Tmin, na.rm = TRUE))
weatherSummary_rakhigarhi$minTemperature.error <- c(
weatherSummary_rakhigarhi$minTemperature.error,
qt(0.975, length(tempData$Tmin) - 1) *
sd(tempData$Tmin, na.rm = TRUE) /
sqrt(length(tempData$Tmin)))
## daily lower and upper deviation
weatherSummary_rakhigarhi$temperature.lowerDeviation <- c(
weatherSummary_rakhigarhi$temperature.lowerDeviation,
mean(tempData$T2M - tempData$Tmin))
weatherSummary_rakhigarhi$temperature.lowerDeviation.error <- c(
weatherSummary_rakhigarhi$temperature.lowerDeviation.error,
qt(0.975, length(tempData$Tmin) - 1) *
sd(tempData$T2M - tempData$Tmin, na.rm = TRUE) /
sqrt(length(tempData$Tmin)))
weatherSummary_rakhigarhi$temperature.upperDeviation <- c(
weatherSummary_rakhigarhi$temperature.upperDeviation,
mean(tempData$Tmax - tempData$T2M))
weatherSummary_rakhigarhi$temperature.upperDeviation.error <- c(
weatherSummary_rakhigarhi$temperature.upperDeviation.error,
qt(0.975, length(tempData$Tmax) - 1) *
sd(tempData$Tmax - tempData$T2M, na.rm = TRUE) /
sqrt(length(tempData$Tmax)))
# precipitation
weatherSummary_rakhigarhi$precipitation.mean <- c(
weatherSummary_rakhigarhi$precipitation.mean,
mean(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.max <- c(
weatherSummary_rakhigarhi$precipitation.max,
max(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.min <- c(
weatherSummary_rakhigarhi$precipitation.min,
min(tempData$RAIN, na.rm = TRUE))
weatherSummary_rakhigarhi$precipitation.error <- c(
weatherSummary_rakhigarhi$precipitation.error,
qt(0.975, length(tempData$RAIN) - 1) *
sd(tempData$RAIN, na.rm = TRUE) /
sqrt(length(tempData$RAIN)))
# ARID
weatherSummary_rakhigarhi$ARID.mean <- c(
weatherSummary_rakhigarhi$ARID.mean,
mean(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.max <- c(
weatherSummary_rakhigarhi$ARID.max,
max(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.min <- c(
weatherSummary_rakhigarhi$ARID.min,
min(tempWatbalData$ARID, na.rm = TRUE))
weatherSummary_rakhigarhi$ARID.error <- c(
weatherSummary_rakhigarhi$ARID.error,
qt(0.975,
length(tempWatbalData$ARID) - 1) *
sd(tempWatbalData$ARID, na.rm = TRUE) /
sqrt(length(tempWatbalData$ARID)))
}
rm(day, tempData, tempWatbalData)
weatherSummary_rakhigarhi <- data.frame(weatherSummary_rakhigarhi)
Get ranges of weather variables:
roundToMultiple <- function(i, baseOfMultiple, roundFunction = round)
{
return(match.fun(roundFunction)(i/baseOfMultiple) * baseOfMultiple)
}
rangeSolar = c(
roundToMultiple(min(c(weatherSummary$solarRadiation.min, weatherSummary_rakhigarhi$solarRadiation.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$solarRadiation.max, weatherSummary_rakhigarhi$solarRadiation.max)), 5, ceiling))
rangeTemp = c(
roundToMultiple(min(c(weatherSummary$minTemperature.min, weatherSummary_rakhigarhi$minTemperature.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$maxTemperature.max, weatherSummary_rakhigarhi$maxTemperature.max)), 5, ceiling))
rangePrecip = c(
roundToMultiple(min(c(weatherSummary$precipitation.min, weatherSummary_rakhigarhi$precipitation.min)), 5, floor),
roundToMultiple(max(c(weatherSummary$precipitation.max, weatherSummary_rakhigarhi$precipitation.max)), 5, ceiling))
Set colours for maximum and minimum temperature:
maxTemperatureColour = hsv(7.3/360, 74.6/100, 70/100)
minTemperatureColour = hsv(232/360, 64.6/100, 73/100)
Define auxiliar function for calculating cumulative curve (extracted from weatherModel.R; see documentation folder in Weather model):
getCumulativePrecipitationOfYear <- function(dailyPrecipitationYear)
{
return(cumsum(dailyPrecipitationYear) / sum(dailyPrecipitationYear))
}
Create figure:
yearLengthInDays_sim = nlevels(factor(weatherSummary$dayOfYear))
yearLengthInDays_real = nlevels(factor(weatherSummary_rakhigarhi$dayOfYear))
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_weatherAndARID.png"
grScale = 2
fontRescale = 0
axisTextRescale = 0
marginTextRescale = 0
leftPlotMargin = 2
rightPlotMargin = 4
cumPrecipMargins <- c(8, 1.5)
bottomMarginWeight = 2
cropRangesHeadLength = 0.1
cropRangesLineWidth = 3
cropRangesLabelsVerticalSpread = 14
cropRangesLabelsSize = 0.95
png(plotName, width = grScale * 800, height = grScale * 1000)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_weatherAndARID.eps"
grScale = 1.2
fontRescale = 0.1
axisTextRescale = -0.1
marginTextRescale = -0.5
leftPlotMargin = 1.5
rightPlotMargin = 3.8
cumPrecipMargins <- c(6, 1.5)
bottomMarginWeight = 1
cropRangesHeadLength = 0.08
cropRangesLineWidth = 3
cropRangesLabelsVerticalSpread = 10
cropRangesLabelsSize = 0.9
extrafont::loadfonts(device = "postscript")
grDevices::cairo_ps(filename = plotName,
pointsize = 12,
width = grScale * 8,
height = grScale * 10,
onefile = FALSE,
fallback_resolution = 600,
family = "sans"
)
}
nColumns = 3
nRowsExceptBottom = 5
layout(rbind(matrix(1:(nColumns * nRowsExceptBottom),
nrow = nRowsExceptBottom,
ncol = nColumns,
byrow = FALSE),
c((nColumns * nRowsExceptBottom) + 1,
rep(((nColumns * nRowsExceptBottom) + 2), 2))),
widths = c(1, 10, 12),
heights = c(rep(10, nRowsExceptBottom) + c(1, rep(0, nRowsExceptBottom - 2), -2), bottomMarginWeight)
)
yLabs <- c(expression(paste("solar radiation (", MJ/m^-2, ")")),
"temperature (C)", "precipitation (mm)", "ARID", "crop seasons")
par(cex = grScale,
cex.axis = grScale * (0.8 + axisTextRescale))
# First column: y axis titles
par(mar = c(0, 0, 0, 0.4))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[1])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.5, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.4, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[3])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.6, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[4])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.6, font = 4,
cex = grScale * (0.78 + fontRescale),
srt = 90,
labels = yLabs[5])
#==============================================================================
# Second column: real-world data
# 1. solar radiation
par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$solarRadiation.mean,
axes = FALSE,
ylim = rangeSolar,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$solarRadiation.mean + weatherSummary_rakhigarhi$solarRadiation.error),
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$solarRadiation.mean - weatherSummary_rakhigarhi$solarRadiation.error),
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$solarRadiation.max,
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$solarRadiation.min,
rev(weatherSummary_rakhigarhi$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangeSolar[1], rangeSolar[2], 5))
# 2. temperature (daily mean, max, min)
## daily mean
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$temperature.mean,
axes = FALSE,
ylim = rangeTemp,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$temperature.mean + weatherSummary_rakhigarhi$temperature.error),
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$temperature.mean - weatherSummary_rakhigarhi$temperature.error),
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$temperature.max,
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$temperature.min,
rev(weatherSummary_rakhigarhi$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
## daily max
lines(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$maxTemperature.mean,
lwd = grScale, col = maxTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$maxTemperature.mean + weatherSummary_rakhigarhi$maxTemperature.error),
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$maxTemperature.mean - weatherSummary_rakhigarhi$maxTemperature.error),
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$maxTemperature.max,
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$maxTemperature.min,
rev(weatherSummary_rakhigarhi$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
## daily min
lines(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$minTemperature.mean,
lwd = grScale, col = minTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$minTemperature.mean + weatherSummary_rakhigarhi$minTemperature.error),
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$minTemperature.mean - weatherSummary_rakhigarhi$minTemperature.error),
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$minTemperature.max,
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$minTemperature.min,
rev(weatherSummary_rakhigarhi$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangeTemp[1], rangeTemp[2], 5))
# 3. precipitation
# cumulative curve
par(mar = c(cumPrecipMargins[1], leftPlotMargin, 0.1, 0.1))
plot(c(1, yearLengthInDays_real), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (year in unique(weatherData_rakhigarhi$year))
{
lines(1:nrow(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year,]),
getCumulativePrecipitationOfYear(weatherData_rakhigarhi[weatherData_rakhigarhi$year == year, "RAIN"]),
lwd = grScale,
col = rgb(0, 0, 0, alpha = 0.05))
}
rm(year)
# daily values
par(new = T,
mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$precipitation.mean,
axes = FALSE,
ylim = rangePrecip * c(1, cumPrecipMargins[2]),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$precipitation.mean + weatherSummary_rakhigarhi$precipitation.error),
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$precipitation.mean - weatherSummary_rakhigarhi$precipitation.error),
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$precipitation.max,
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$precipitation.min,
rev(weatherSummary_rakhigarhi$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(rangePrecip[1], rangePrecip[2], 10))
# 4. ARID
par(mar = c(0.1, leftPlotMargin, 0.1, 0.1))
plot(1:yearLengthInDays_real,
weatherSummary_rakhigarhi$ARID.mean,
axes = FALSE,
ylim = c(0, 1),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$ARID.mean + weatherSummary_rakhigarhi$ARID.error),
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c((weatherSummary_rakhigarhi$ARID.mean - weatherSummary_rakhigarhi$ARID.error),
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$ARID.max,
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_real,
rev(1:yearLengthInDays_real)),
y = c(weatherSummary_rakhigarhi$ARID.min,
rev(weatherSummary_rakhigarhi$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
axis(2, at = seq(0, 1, 0.1))
#5. crop growing seasons
par(mar = c(3, leftPlotMargin, 0.1, 0.1))
plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (aCrop in levels(cropTable$crop))
{
tempData <- cropTable[cropTable$crop == aCrop,]
aCropIndex <- match(aCrop, levels(cropTable$crop))
# winter crop
if (tempData$sugSowingDay > tempData$sugHarvestingDay)
{
arrows(x0 = 1, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = 365, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
else
{
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
}
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# x axis
axis(1, at = c(0,
31,
31+28,
31+28+31,
31+28+31+30,
31+28+31+30+31,
31+28+31+30+31+30,
31+28+31+30+31+30+31,
31+28+31+30+31+30+31+31,
31+28+31+30+31+30+31+31+30,
31+28+31+30+31+30+31+31+30+31,
31+28+31+30+31+30+31+31+30+31+30,
31+28+31+30+31+30+31+31+30+31+30+31),
las = 2
)
#==============================================================================
# Third column: simulated data
# 1. solar radiation
par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$solarRadiation.mean,
axes = FALSE,
ylim = rangeSolar,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$solarRadiation.mean + weatherSummary$solarRadiation.error),
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$solarRadiation.mean - weatherSummary$solarRadiation.error),
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$solarRadiation.max,
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$solarRadiation.min,
rev(weatherSummary$solarRadiation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 2. temperature (daily mean, max, min)
## daily mean
plot(1:yearLengthInDays_sim,
weatherSummary$temperature.mean,
axes = FALSE,
ylim = rangeTemp,
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$temperature.mean + weatherSummary$temperature.error),
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$temperature.mean - weatherSummary$temperature.error),
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$temperature.max,
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$temperature.min,
rev(weatherSummary$temperature.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
## daily max
lines(1:yearLengthInDays_sim,
weatherSummary$maxTemperature.mean,
lwd = grScale, col = maxTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$maxTemperature.mean + weatherSummary$maxTemperature.error),
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$maxTemperature.mean - weatherSummary$maxTemperature.error),
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$maxTemperature.max,
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$maxTemperature.min,
rev(weatherSummary$maxTemperature.mean)),
col = adjustcolor(maxTemperatureColour, alpha.f = 0.3),
border = NA)
## daily min
lines(1:yearLengthInDays_sim,
weatherSummary$minTemperature.mean,
lwd = grScale, col = minTemperatureColour)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$minTemperature.mean + weatherSummary$minTemperature.error),
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$minTemperature.mean - weatherSummary$minTemperature.error),
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$minTemperature.max,
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$minTemperature.min,
rev(weatherSummary$minTemperature.mean)),
col = adjustcolor(minTemperatureColour, alpha.f = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 3. precipitation
# cumulative curve
par(mar = c(cumPrecipMargins[1], 0.1, 0.1, rightPlotMargin))
plot(c(1, yearLengthInDays_sim), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (randomSeed in unique(weatherData$randomSeed))
{
for (year in unique(weatherData$currentYear))
{
lines(1:nrow(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year,]),
getCumulativePrecipitationOfYear(weatherData[weatherData$randomSeed == randomSeed & weatherData$currentYear == year, "RAIN"]),
lwd = grScale,
col = rgb(0, 0, 0, alpha = 0.05))
}
}
rm(year, randomSeed)
axis(4, at = seq(0, 1, 0.25))
mtext("cumulative annual sum", 4, line = 2,
cex = grScale * (1.5 + marginTextRescale))
# daily values
par(new = T,
mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$precipitation.mean,
axes = FALSE,
ylim = rangePrecip * c(1, cumPrecipMargins[2]),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$precipitation.mean + weatherSummary$precipitation.error),
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$precipitation.mean - weatherSummary$precipitation.error),
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$precipitation.max,
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$precipitation.min,
rev(weatherSummary$precipitation.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# 4. ARID
par(mar = c(0.1, 0.1, 0.1, rightPlotMargin))
plot(1:yearLengthInDays_sim,
weatherSummary$ARID.mean,
axes = FALSE,
ylim = c(0, 1),
type = "l", lwd = grScale)
## 95% confidence interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$ARID.mean + weatherSummary$ARID.error),
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c((weatherSummary$ARID.mean - weatherSummary$ARID.error),
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.5),
border = NA)
## min-max interval
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$ARID.max,
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border=NA)
polygon(x = c(1:yearLengthInDays_sim,
rev(1:yearLengthInDays_sim)),
y = c(weatherSummary$ARID.min,
rev(weatherSummary$ARID.mean)),
col = rgb(0,0,0, alpha = 0.3),
border = NA)
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
#5. crop growing seasons
par(mar = c(3, 0.1, 0.1, rightPlotMargin))
plot(c(1, 365), c(1, nrow(cropTable)), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
for (aCrop in levels(cropTable$crop))
{
tempData <- cropTable[cropTable$crop == aCrop,]
aCropIndex <- match(aCrop, levels(cropTable$crop))
# winter crop
if (tempData$sugSowingDay > tempData$sugHarvestingDay)
{
arrows(x0 = 1, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = yearLengthInDays_sim, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
else
{
arrows(x0 = tempData$sugSowingDay, y0 = aCropIndex,
x1 = tempData$sugHarvestingDay, y1 = aCropIndex,
angle = 90, code = 3, col = cropColours[aCropIndex],
length = cropRangesHeadLength, lwd = cropRangesLineWidth * grScale)
}
mtext(aCrop,
side = 4, col = cropColours[aCropIndex],
padj = ((1 - aCropIndex/nrow(cropTable)) - 0.4) * cropRangesLabelsVerticalSpread,
las = 2, cex = grScale * cropRangesLabelsSize)
}
# solstices
abline(v = 31+28+31+30+31+21, # 21 June (approx.)
lty = 3, lwd = grScale)
abline(v = 31+28+31+30+31+30+31+31+30+31+30+21, # 21 December (approx.)
lty = 3, lwd = grScale)
# x axis
axis(1, at = c(0,
31,
31+28,
31+28+31,
31+28+31+30,
31+28+31+30+31,
31+28+31+30+31+30,
31+28+31+30+31+30+31,
31+28+31+30+31+30+31+31,
31+28+31+30+31+30+31+31+30,
31+28+31+30+31+30+31+31+30+31,
31+28+31+30+31+30+31+31+30+31+30,
31+28+31+30+31+30+31+31+30+31+30+31),
las = 2
)
# bottom row: empty and "day of year" or x axis title
par(mar = c(0, 0, 0, 0))
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.45, y = 0.5, font = 4,
cex = grScale * (0.8 + fontRescale),
labels = "day of year")
dev.off()
}
rm(imageFormat, grScale, fontRescale, yLabs, nColumns, nRowsExceptBottom,
axisTextRescale, marginTextRescale, leftPlotMargin, rightPlotMargin, cumPrecipMargins, bottomMarginWeight,
tempData, aCrop, aCropIndex, cropRangesHeadLength, cropRangesLineWidth)
knitr::include_graphics(plotName)
Display weather parameters:
knitr::kable(t(yieldData[1, 2:26]), col.names = "default (NW India echo)")
| default (NW India echo) | |
|---|---|
| temperature_annualMaxAt2m | 37.00 |
| temperature_annualMinAt2m | 12.80 |
| temperature_meanDailyFluctuation | 2.20 |
| temperature_dailyLowerDeviation | 6.80 |
| temperature_dailyUpperDeviation | 7.90 |
| solar_annualMax | 24.20 |
| solar_annualMin | 9.20 |
| solar_meanDailyFluctuation | 3.30 |
| CO2_annualMin | 245.00 |
| CO2_annualMax | 255.00 |
| CO2_meanDailyFluctuation | 1.00 |
| precipitation_yearlyMean | 489.00 |
| precipitation_yearlySd | 142.20 |
| precipitation_dailyCum_nSamples | 200.00 |
| precipitation_dailyCum_maxSampleSize | 10.00 |
| precipitation_dailyCum_plateauValue_yearlyMean | 0.25 |
| precipitation_dailyCum_plateauValue_yearlySd | 0.10 |
| precipitation_dailyCum_inflection1_yearlyMean | 40.00 |
| precipitation_dailyCum_inflection1_yearlySd | 5.00 |
| precipitation_dailyCum_rate1_yearlyMean | 0.07 |
| precipitation_dailyCum_rate1_yearlySd | 0.02 |
| precipitation_dailyCum_inflection2_yearlyMean | 240.00 |
| precipitation_dailyCum_inflection2_yearlySd | 20.00 |
| precipitation_dailyCum_rate2_yearlyMean | 0.08 |
| precipitation_dailyCum_rate2_yearlySd | 0.02 |
Summary statistics per crop:
yieldData_summary <- tapply(as.numeric(as.character(yieldData$yield)), yieldData$crop, summary)
yieldData_summary <- data.frame(Reduce(rbind, yieldData_summary), row.names = names(yieldData_summary))
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
## Warning in f(init, x[[i]]): number of columns of result is not a multiple of
## vector length (arg 2)
knitr::kable(yieldData_summary)
| Min. | X1st.Qu. | Median | Mean | X3rd.Qu. | Max. | |
|---|---|---|---|---|---|---|
| proso millet | 885.6262 | 1059.3122 | 1109.8132 | 1108.5131 | 1157.0446 | 1373.1612 |
| pearl millet | 300.5779 | 332.6118 | 342.6917 | 342.9984 | 352.8890 | 392.0383 |
| rice | 0.0000 | 230.7686 | 265.1131 | 261.6839 | 299.6243 | 397.6978 |
| barley | 337.8926 | 389.9014 | 412.2040 | 414.7265 | 437.7929 | 520.7963 |
| wheat 1 | 280.4535 | 327.1036 | 345.8900 | 347.8794 | 367.3557 | 437.3639 |
| wheat 2 | 296.1396 | 338.7688 | 357.1001 | 359.2555 | 378.5212 | 447.0028 |
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldPerCrop.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldPerCrop.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
par(mar = c(6,6,1,1), cex.lab = grScale, cex.axis = 0.8 * grScale)
boxplot(yield ~ factor(crop), data = yieldData,
ylab = expression(paste("yield (", g/m^2, ")")), xlab = "",
col = cropColours
)
dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)
ggplot2 option:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldPerCrop_ggplot.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldPerCrop_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
print({
ggplot(yieldData, aes(y=yield, x=crop, fill=crop)) +
geom_violin(position="dodge") +
#coord_flip() +
#geom_boxplot() +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
scale_x_discrete() +
theme_light(base_size = 11 * grScale) +
ylab(expression(paste("yield (", g/m^2, ")")))
})
dev.off()
}
## Warning: Removed 75 rows containing non-finite values (stat_ydensity).
## Warning: Removed 75 rows containing non-finite values (stat_ydensity).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow.png"
grScale = 2
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
layout(matrix(c(1, 2), nrow = 1, ncol = 2),
widths = c(10, 3))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
boxplot(meanARID_grow ~ factor(crop), data = yieldData,#[yieldData$yield > 0,], # show only non-zero yield
ylab = "mean ARID", xlab = "growing seasons",
#las = 2,
col = cropColours)
text(x = 0.5, y = minARID * 1.005, labels = "a", font = 2, cex = grScale)
par(mar = c(5,1,1,1))
xhist <- hist(yieldData$meanARID, breaks = 20, plot = FALSE)
barplot(xhist$counts, axes = FALSE, space = 0, horiz = TRUE,
xlab = "current year")
text(x = max(xhist$counts) * 0.9, y = minARID * 1.005, labels = "b", font = 2, cex = grScale)
dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)
ggplot2 option:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow_ggplot.png"
grScale = 2
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_ARIDVsARID_grow_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 6 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
tempData <- rbind(yieldData[, c("crop", "meanARID_grow")],
cbind(crop = rep("year", nrow(yieldData)),
meanARID_grow = yieldData$meanARID))
tempData$meanARID_grow <- as.numeric(tempData$meanARID_grow)
print({
(ggplot(tempData, aes(fill = crop, y = meanARID_grow, x = crop)) +
geom_violin(position = "dodge") +
#geom_boxplot() +
scale_fill_manual(values = c(cropColours, "grey"), name="", guide = FALSE) +
scale_x_discrete() +
theme_light(base_size = 11 * grScale) +
theme(axis.title.x = element_blank()) +
ylab("mean ARID")) +
ylim(minARID, maxARID)
})
dev.off()
}
## Warning: Removed 76 rows containing non-finite values (stat_ydensity).
## Warning: Removed 76 rows containing non-finite values (stat_ydensity).
rm(tempData, imageFormat, grScale)
knitr::include_graphics(plotName)
Plotting absolute values of yield vs ARID:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldVsARID.png"
grScale = 2
cex.points = 1
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldVsARID.eps"
grScale = 2
cex.points = 0.5
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
widths = c(10, 3.5))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
plot(c(minARID, maxARID),
c(minYield, maxYield),
xlab = "mean ARID (harvest year)",
ylab = expression(paste("yield (", g/m^2, ")")))
for (aCrop in levels(yieldData$crop))
{
points(yieldData$meanARID[yieldData$crop == aCrop],
yieldData$yield[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(yield ~ meanARID, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minARID * 1.005, y = maxYield * 0.95, labels = "a", font = 2, cex = grScale)
plot(c(minARID, maxARID),
c(minYield, maxYield),
xlab = "mean ARID (grow season)",
ylab = expression(paste("yield (", g/m^2, ")")))
for (aCrop in levels(yieldData$crop))
{
points(yieldData$meanARID_grow[yieldData$crop == aCrop],
yieldData$yield[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minARID * 1.005, y = maxYield * 0.95, labels = "b", font = 2, cex = grScale)
par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0, y = 1,
legend = stringi::stri_c(cropTable$S_water[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"),
title = "S_water",
fill = cropColours)
dev.off()
}
rm(imageFormat, grScale)
knitr::include_graphics(plotName)
ggplot alternative:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp1_yieldVsARID_ggplot.png"
grScale = 2
png(plotName, width = grScale * 350, height = grScale * 500)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp1_yieldVsARID_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 3.5 * grScale,
height = 5 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
print({
ggplot(data = yieldData, aes(x = meanARID_grow, y = yield, color = crop, fill = crop)) +
geom_point(size = 1) +
geom_smooth(method = 'lm') +
scale_color_manual(values = cropColours, name="", guide = FALSE) +
scale_fill_manual(values = cropColours, name="", guide = FALSE) +
facet_wrap( ~ crop, ncol = 2, scales = 'free') +
xlab("mean ARID") +
ylab(expression(paste("yield (", g/m^2, ")"))) +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))
})
dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)
Calculate correlation and linear models fitness:
linearModels.table <- data.frame(
levels(yieldData$crop),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop))
)
names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(yieldData$crop))
{
linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <-
cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop], y = yieldData$yield[yieldData$crop == aCrop], use = "pairwise.complete.obs")
linearModel <- lm(yield ~ meanARID_grow, data = yieldData[yieldData$crop == aCrop,])
linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between ARID (x) and Yield (y) per crop")
| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.0201669 | 1086.7516 | 0 | 23.79756 | 0.5813393 | -0.0009296 |
| pearl millet | -0.0482022 | 347.6334 | 0 | -7.44716 | 0.1872903 | 0.0009897 |
| rice | -0.9094099 | 693.8981 | 0 | -637.77881 | 0.0000000 | 0.8267950 |
| barley | -0.8629164 | 629.6293 | 0 | -273.61862 | 0.0000000 | 0.7442715 |
| wheat 1 | -0.8470997 | 531.3303 | 0 | -231.94635 | 0.0000000 | 0.7171872 |
| wheat 2 | -0.8582940 | 540.0265 | 0 | -228.55794 | 0.0000000 | 0.7363044 |
Clear temporary objects:
rm(linearModel, linearModels.table,
watbalOutput_rakhigarhi, watbalParameters, soil.FC, soil.WP,
weatherData, weatherData_rakhigarhi, weatherSummary, weatherSummary_rakhigarhi, weatherSummaryNames,
xhist, yieldData, yieldData_summary,
yearLengthInDays_real, yearLengthInDays_sim,
maxARID, minARID, maxYield, minYield, rangePrecip, rangeSolar, rangeTemp,
maxTemperatureColour, minTemperatureColour,
plotName)
## Warning in rm(linearModel, linearModels.table, watbalOutput_rakhigarhi, : object
## 'soil.FC' not found
## Warning in rm(linearModel, linearModels.table, watbalOutput_rakhigarhi, : object
## 'soil.WP' not found
Experiment 2 perform the same runs of experiment 1, except by varying stochastically two weather model parameters concerning precipitation: precipitation_yearlyMean (average annual sum of precipitation) and precipitation_dailyCum_plateauValue_yearlyMean (average plateau value of the cumulative curve of daily precipitation).
Load data:
yieldData <- read.csv("models/output/SIMPLE-crop-model_yield-exp_type-of-experiment=precipitation-variation_experiment-name=exp2_initRandomSeed=0.csv")
Convert yield = 0 to NA when meanARID_grow = NA (i.e., first year simulation runs with incomplete growing seasons):
yieldData$yield[is.na(yieldData$meanARID_grow)] <- NA
Force the same order of crops used by cropTable:
yieldData$crop <- factor(yieldData$crop, levels = levels(cropTable$crop))
Calculate the range of ARID, yield, precipitation annual sum (yearly mean) and the “plateau value” (yearly mean) of the cumulative curve of daily precipitation:
minARID = round(min(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
maxARID = round(max(c(yieldData$meanARID, yieldData$meanARID_grow), na.rm = TRUE), digits = 2)
minYield = round(min(yieldData$yield, na.rm = TRUE), digits = -1)
maxYield = round(max(yieldData$yield, na.rm = TRUE), digits = -1)
minPrecipitation_yearlyMean = round(min(yieldData$precipitation_yearlyMean), digits = -2)
maxPrecipitation_yearlyMean = round(max(yieldData$precipitation_yearlyMean), digits = -2)
minPrecipitation_dailyCum_plateauValue_yearlyMean = round(min(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)
maxPrecipitation_dailyCum_plateauValue_yearlyMean = round(max(yieldData$precipitation_dailyCum_plateauValue_yearlyMean), digits = 2)
Display explored values of precipitation_yearlyMean and precipitation_dailyCum_plateauValue_yearlyMean:
layout(matrix(1:2, ncol = 2))
hist(yieldData$precipitation_yearlyMean, main = "precipitation_yearlyMean")
hist(yieldData$precipitation_dailyCum_plateauValue_yearlyMean, main = "precipitation_dailyCum_plateauValue_yearlyMean")
Plotting ARID vs precipitation yearly mean and seasonal balance (plateau value):
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation.png"
grScale = 2
cex.points = 1
png(plotName, width = grScale * 500, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation.eps"
grScale = 2
cex.points = 0.5
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
layout(matrix(c(1, 2, 3, 3), nrow = 2, ncol = 2, byrow = FALSE),
widths = c(10, 3.5))
par(mar = c(5,5,1,1), cex.lab = 0.8 * grScale)
plot(c(minPrecipitation_yearlyMean, maxPrecipitation_yearlyMean),
c(minARID, maxARID),
xlab = "precipitation annual sum (mm) (yearly mean)",
ylab = "mean ARID (growing seasons)")
for (aCrop in levels(yieldData$crop))
{
points(yieldData$precipitation_yearlyMean[yieldData$crop == aCrop],
yieldData$meanARID_grow[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minPrecipitation_yearlyMean * 0.95, y = minARID * 1.2, labels = "a", font = 2, cex = grScale)
plot(c(minPrecipitation_dailyCum_plateauValue_yearlyMean, maxPrecipitation_dailyCum_plateauValue_yearlyMean),
c(minARID, maxARID),
xlab = "cumulative daily precipitation plateau value (yearly mean)",
ylab = "mean ARID (growing seasons)")
for (aCrop in levels(yieldData$crop))
{
points(yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop],
yieldData$meanARID_grow[yieldData$crop == aCrop],
col = cropColours[match(aCrop, levels(yieldData$crop))],
pch = 20, cex = cex.points
)
abline(lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,]),
col = cropColours[match(aCrop, levels(yieldData$crop))],
lwd = 2)
}
rm(aCrop)
text(x = minPrecipitation_dailyCum_plateauValue_yearlyMean * 0.95, y = minARID * 1.2, labels = "b", font = 2, cex = grScale)
par(mar = c(0, 0, 0, 0), cex = 0.8 * grScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0, y = 1,
legend = stringi::stri_c(cropTable$sugSowingDay[order(match(cropTable$crop, levels(cropTable$crop)))], " (", levels(yieldData$crop), ")"),
title = "sowing day of year",
fill = cropColours)
dev.off()
}
rm(imageFormat, grScale, cex.points)
knitr::include_graphics(plotName)
ggplot alternative:
for (imageFormat in listOfImageFormats)
{
if (imageFormat == "png")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation_ggplot.png"
grScale = 2
png(plotName, width = grScale * 600, height = grScale * 300)
}
if (imageFormat == "eps")
{
plotName = "figures/Fig_exp2_ARIDvsPrecipitation_ggplot.eps"
grScale = 2
extrafont::loadfonts(device = "postscript")
grDevices::postscript(file = plotName,
pointsize = 10,
width = 5 * grScale,
height = 3 * grScale,
horizontal = FALSE,
paper = "special",
onefile = FALSE,
family = "sans",
colormodel = "cmyk")
}
print({
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_yearlyMean, color = crop)) +
geom_point(size = 1) +
geom_smooth(method = 'lm') +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("mean ARID") +
xlab("precipitation annual sum (mm) (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
(ggplot(data = yieldData, aes(y = meanARID_grow, x = precipitation_dailyCum_plateauValue_yearlyMean, color = crop)) +
geom_point(size = 1) +
geom_smooth(method = 'lm') +
scale_color_discrete_qualitative(guide = FALSE) +
facet_wrap( ~ crop, ncol = 2) +
ylab("") +
xlab("cumulative daily precipitation plateau value (yearly mean)") +
theme_light(base_size = 11 * grScale) +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 15),
strip.background = element_rect(fill = "white"),
strip.text = element_text(colour = "black"))) +
plot_layout(nrow = 1, widths = c(1, 1))
})
dev.off()
}
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 75 rows containing non-finite values (stat_smooth).
## Warning: Removed 75 rows containing missing values (geom_point).
rm(imageFormat, grScale)
knitr::include_graphics(plotName)
Calculate correlation and linear models fitness:
linearModels.table <- data.frame(
levels(yieldData$crop),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop))
)
names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(yieldData$crop))
{
linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <-
cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop],
y = yieldData$precipitation_yearlyMean[yieldData$crop == aCrop],
use = "pairwise.complete.obs")
linearModel <- lm(meanARID_grow ~ precipitation_yearlyMean, data = yieldData[yieldData$crop == aCrop,])
linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between precipitation annual sum yearly mean (x) and ARID (y) per crop")
| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | -0.5760686 | 0.9960945 | 0 | -0.0001314 | 0 | 0.3309618 |
| pearl millet | -0.7070320 | 0.9428775 | 0 | -0.0003964 | 0 | 0.4992257 |
| rice | -0.7308760 | 0.9545028 | 0 | -0.0003454 | 0 | 0.5335569 |
| barley | -0.7133386 | 0.8883206 | 0 | -0.0004509 | 0 | 0.5081726 |
| wheat 1 | -0.7138806 | 0.8916387 | 0 | -0.0004425 | 0 | 0.5089473 |
| wheat 2 | -0.7138806 | 0.8916387 | 0 | -0.0004425 | 0 | 0.5089473 |
linearModels.table <- data.frame(
levels(yieldData$crop),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop)),
numeric(nlevels(yieldData$crop))
)
names(linearModels.table) <- c("crop", "correlation_Pearson", "intercept", "intercept_p", "speed", "speed_p", "adjrsquared")
for (aCrop in levels(yieldData$crop))
{
linearModels.table$correlation_Pearson[linearModels.table$crop == aCrop] <-
cor(x = yieldData$meanARID_grow[yieldData$crop == aCrop],
y = yieldData$precipitation_dailyCum_plateauValue_yearlyMean[yieldData$crop == aCrop],
use = "pairwise.complete.obs")
linearModel <- lm(meanARID_grow ~ precipitation_dailyCum_plateauValue_yearlyMean, data = yieldData[yieldData$crop == aCrop,])
linearModels.table$intercept[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,1]
linearModels.table$intercept_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[1,4]
linearModels.table$speed[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,1]
linearModels.table$speed_p[linearModels.table$crop == aCrop] <- summary(linearModel)$coefficients[2,4]
linearModels.table$adjrsquared[linearModels.table$crop == aCrop] <- summary(linearModel)$adj.r.squared
}
rm(aCrop)
knitr::kable(linearModels.table, caption = "Pearson's correlation and linear model between cumulative daily precipitation plateau value yearly mean (x) and ARID (y) per crop")
| crop | correlation_Pearson | intercept | intercept_p | speed | speed_p | adjrsquared |
|---|---|---|---|---|---|---|
| proso millet | 0.1633846 | 0.8971123 | 0 | 0.0628231 | 6.9e-06 | 0.0253933 |
| pearl millet | 0.4151101 | 0.5367639 | 0 | 0.3922293 | 0.0e+00 | 0.1712098 |
| rice | 0.4268577 | 0.6015239 | 0 | 0.3400367 | 0.0e+00 | 0.1811142 |
| barley | -0.3736394 | 0.8743004 | 0 | -0.3980747 | 0.0e+00 | 0.1384163 |
| wheat 1 | -0.3754886 | 0.8787469 | 0 | -0.3923104 | 0.0e+00 | 0.1398036 |
| wheat 2 | -0.3754886 | 0.8787469 | 0 | -0.3923104 | 0.0e+00 | 0.1398036 |
Clear temporary objects:
rm(linearModel, linearModels.table,
yieldData,
maxARID, minARID, maxYield, minYield,
maxPrecipitation_yearlyMean, minPrecipitation_yearlyMean,
maxPrecipitation_dailyCum_plateauValue_yearlyMean, minPrecipitation_dailyCum_plateauValue_yearlyMean,
plotName)
Defining terrain random number generator seed (must be available inside “terrains” folder:
listOfTerrainSeeds <- c("0", "35", "56", "72", "92")
Setting dimensions:
figScale = 2
terrainDim = 450 * figScale
topHeight = 50 * figScale
bottomHeight = 150 * figScale
leftWidth = 100 * figScale
internalMargin = 5 * figScale
figWidth = terrainDim * 3 + leftWidth + 2 * internalMargin
figHeight = terrainDim * length(listOfTerrainSeeds) + topHeight + bottomHeight + (length(listOfTerrainSeeds) - 1) * internalMargin
Create margin header images:
headers <- c("seed",
"Elevation",
"Soil texture",
"Soil texture type")
namesInFile <- c("seed", "elev", "soilText", "soilTextType")
for (i in 1:length(headers))
{
#if (i == 1) { width = leftWidth } else { width = terrainDim }
width = ifelse(i == 1, leftWidth, terrainDim)
png(paste0("FigTerrains_tempFiles/FigTerrains_", namesInFile[i], ".PNG"),
width = width,
height = topHeight)
par(mar = rep(0, 4))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.7, cex = figScale * 3.7, adj = c(0.5, 0.5),
labels = headers[i])
dev.off()
}
rm(i)
for (seed in listOfTerrainSeeds)
{
png(paste0("FigTerrains_tempFiles/FigTerrains_seed=", seed, ".PNG"),
width = leftWidth,
height = terrainDim)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.5, cex = figScale * 5, adj = .5,
labels = seed)
dev.off()
}
rm(seed)
Create legend images:
# code converted from NetLogo (maximum = 255)
elev <- 100/255 + (155/255) * seq(0, 1, length.out = 10)
png("FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = c(2,2,1,2), cex = figScale * 2)
barplot(rep(100, 10),
col = rgb(elev - (100/255), elev, 0),
border = rgb(elev - (100/255), elev, 0),
axes = FALSE, space = 0)
mtext("min.", side = 1, line = 0.1, adj = 0, cex = figScale * 2)
mtext("max.", side = 1, line = 0.1, adj = 1, cex = figScale * 2)
#mtext("elevation", side = 1, line = 1.2, adj = .5, cex = figScale * 3)
dev.off()
## png
## 2
#======================================================================
soilTexture <- c("100% sand", "100% silt", "100% clay")
soilTexture_colors <- c("red", "green", "blue")
xy <- list(c(0.2, 1.05),
c(0.2, 0.75),
c(0.2, 0.45))
png("FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTexture))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTexture[i],
fill = soilTexture_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 2.5, pt.cex = figScale * 5)
}
rm(i)
dev.off()
## png
## 2
#======================================================================
soilTextureTypes <- c("Sand", "Loamy sand", "Sandy loam",
"Loam", "Silt loam", "Silty clay",
"Clay", "Clay loam", "Sandy clay loam")
soilTextureTypes_colors <- c("#d73229", "#f16a15", "#9d6e48",
"#eded31", "#59b03c", "#54c4c4",
"#2d8dbe", "#345da9", "#a71b6a")
xy <- list(c(-0.08, 1), c(0.17, 1), c(0.55, 1),
c(-0.08, 0.7), c(0.17, 0.7), c(0.55, 0.7),
c(-0.08, 0.4), c(0.17, 0.4), c(0.55, 0.4))
png("FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTextureTypes))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTextureTypes[i],
fill = soilTextureTypes_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 1.8, pt.cex = figScale * 5)
}
rm(i)
dev.off()
## png
## 2
Create blank image:
png("FigTerrains_tempFiles/FigTerrains_blank.PNG",
width = 10,
height = 10)
## Warning in png("FigTerrains_tempFiles/FigTerrains_blank.PNG", width = 10, :
## 'width=10, height=10' are unlikely values in pixels
dev.off()
## png
## 2
Load input images:
elevTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,".PNG")
soilTextTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_soilTexture.PNG")
soilTextTypeTerrains <- paste0("models/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfTerrainSeeds ,"_soilTextureTypes.PNG")
numberOfColumns = 6
numberOfRows = 11
# it is important that these plots were saved as "PNG not "png"
listOfPlots <- c("FigTerrains_tempFiles/FigTerrains_seed.PNG",
"FigTerrains_tempFiles/FigTerrains_elev.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_soilText.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_soilTextType.PNG",
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[1], ".PNG"),
elevTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[1],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[1],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[2], ".PNG"),
elevTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[2],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[2],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[3], ".PNG"),
elevTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[3],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[3],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[4], ".PNG"),
elevTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[4],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[4],
# ---
rep("FigTerrains_tempFiles/FigTerrains_blank.PNG", numberOfColumns),
# ---
paste0("FigTerrains_tempFiles/FigTerrains_seed=", listOfTerrainSeeds[5], ".PNG"),
elevTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTerrains[5],
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
soilTextTypeTerrains[5],
# ---
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_elev.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_soilText.PNG",
"FigTerrains_tempFiles/FigTerrains_blank.PNG",
"FigTerrains_tempFiles/FigTerrains_legend_soilTextType.PNG"
)
grobs <- lapply(lapply(listOfPlots, png::readPNG), grid::rasterGrob)
Create Figure:
widths = 50 * c(leftWidth/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth) * figScale
heights = (figHeight/figWidth) * 50 * c(topHeight/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
bottomHeight/figHeight) * figScale
lay <- matrix(1:(numberOfColumns * numberOfRows), ncol = numberOfColumns, byrow = TRUE)
png("figures/Fig_terrains.png", width = figWidth, height = figHeight)
gridExtra::grid.arrange(grobs = grobs,
layout_matrix = lay,
widths = unit(widths, "cm"),
heights = unit(heights, "cm"))
dev.off()
## png
## 2
knitr::include_graphics("figures/Fig_terrains.png")